%--------------------------------------------------------------------------
% 脚本名称：pipeline1_3DPSF
% 用途：采用噪声原图得到潜在清晰图像，在某一个聚焦波段位置上重建PSF
% 输出：在输出文件夹中输出PSF图像，如果有151个波段，则有151张图像。且有一个
%       mat文件保存各波段的模糊核
% 作者：刘昀卓 522543417@qq.com
% 时间：2022/3/30
%--------------------------------------------------------------------------

clc;clear;close;
% 添加必要库
addpath(genpath('../utils'));

%-----------------------------指定文件路径及参数---------------------------
SAVE_FOLDER = '结果'; % 结果保存路径
INPUT_FOLDER = 'F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据'; % 输入文件夹路径
FOCUSSHIFT_WAVELENGTH = 550; %nm 重建的三维PSF是哪个波长的聚焦位置
FOCUS_RADIUS = 30; %像素 聚焦波段圆点图像半径
DEFOCUS_WAVELENGTH = 650; %nm 模糊核大小估计所用的离焦波段
DEFOCUS_RADIUS = 48; %像素 离焦波段圆点图像半径
MAX_PSF_SIZE = 121; % 最大的模糊核估计尺寸 也是填充显示的PSF图像的尺寸
WAVELENGTH_RANGE = 420:10:720; % 要重建的波长范围
EXCEL_PATH = 'F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\psf数据\版画采光谱图像位置及曝光时间数据.xlsx';% 输入的各波段聚焦位置EXCEL表
% RECTANGLE_ROI = [670,200,720,540]; % 用于裁剪图像的矩形框
savearray_path = ""; % 为空代表重新从新开始重建
cb_rows = 15;
cb_cols = 20;
rect = [278,187,174,125];
% 原图
cb_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\cb.png"; % 棋盘格原图
noise_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\noise1.png"; % 噪声原图
black_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\black.png"; % 黑色原图
white_prime = "F:\OneDrive\lyz资料\程序资料\软件相关\前处理\Step5 PSF重建\原始图像淘宝\white.png"; % 白色原图
%--------------------------------------------------------------------------

%% 1.圆点目标物拍摄图测量其半径(相似原理确定核R大小)
% 1.1 聚焦搜索半径
focus_radius = FOCUS_RADIUS;
% 1.2 离焦搜索半径
defocus_radius = DEFOCUS_RADIUS;
% 1.3 通过相似原理得到比例放大系数 来为后续核估计大小
ExcelFiledata = xlsread(EXCEL_PATH,'Sheet1');% 从Excel表中导出所有数值数据
All_Wavelength = ExcelFiledata(1:end,1)';% 取出All_Wavelength,拍摄波段
All_Position = ExcelFiledata(1:end,2)';% 取出All_Position,聚焦位置
focus_shift =  All_Position(All_Wavelength==550); % 550nm聚焦位置
defocus_shift =  All_Position(All_Wavelength==650); % 650nm离焦位置
magnify = abs(focus_radius-defocus_radius)/abs(focus_shift-defocus_shift); 

%% 2.读入伯努利噪声、白板、黑板、棋盘格原图
cb_prime = im2double(imread(cb_prime));
noise_prime = im2double(imread(noise_prime)); % 拍摄的噪声图像
black_prime = im2double(imread(black_prime)); % 拍摄的黑板
white_prime = im2double(imread(white_prime));% 拍摄的白板

% % 裁剪拍摄的图像
% rectangleROI = RECTANGLE_ROI;
% rectangleROI(3) = rectangleROI(3)-1;
% rectangleROI(4) = rectangleROI(4)-1;

%% 3 潜在清晰图像配准
% 3.1聚焦位置数据
filename = fullfile(INPUT_FOLDER,'cb','光谱图像',[num2str(FOCUSSHIFT_WAVELENGTH),'nm(',num2str(All_Position(All_Wavelength==550)),'mm).tif']);
img_cb = im2double(imread(filename));
% 3.2 进行图像配准
iprime_norm = warp_function(cb_prime,...
                noise_prime,...
                img_cb,...
                cb_rows,...
                cb_cols);

%% 4 逐波段重建PSF
savepath = SAVE_FOLDER;
imgpath = fullfile(savepath,'PSF图片');
img1path = fullfile(savepath,'去模糊图片');
if ~(exist(savepath,'dir')==7)
   mkdir(savepath);
end
if ~(exist(imgpath,'dir')==7)
   mkdir(imgpath);
end
if ~(exist(img1path,'dir')==7)
   mkdir(img1path);
end

psf_size = MAX_PSF_SIZE;
% 加载保存的savearray
if savearray_path == ""
    savearray = cell(length(WAVELENGTH_RANGE),1);
else
    load(savearray_path)
end

for i=1:length(WAVELENGTH_RANGE)
    % 4.1 读入拍摄的伯努利噪声、白板、黑板、棋盘格
    % 离焦位置数据
    wavelength = WAVELENGTH_RANGE(i);
    img_noise = imread(fullfile(INPUT_FOLDER, 'noise1' ,[num2str(FOCUSSHIFT_WAVELENGTH),'nm聚集位置所有中心波长图像'],[num2str(wavelength),'nm(',num2str(All_Position(All_Wavelength==FOCUSSHIFT_WAVELENGTH)),'mm).tif']));
    img_noise = im2double(img_noise);
    img_blk = imread(fullfile(INPUT_FOLDER, 'black' ,[num2str(FOCUSSHIFT_WAVELENGTH),'nm聚集位置所有中心波长图像'],[num2str(wavelength),'nm(',num2str(All_Position(All_Wavelength==FOCUSSHIFT_WAVELENGTH)),'mm).tif']));
    img_blk = im2double(img_blk);
    img_wht = imread(fullfile(INPUT_FOLDER, 'white' ,[num2str(FOCUSSHIFT_WAVELENGTH),'nm聚集位置所有中心波长图像'],[num2str(wavelength),'nm(',num2str(All_Position(All_Wavelength==FOCUSSHIFT_WAVELENGTH)),'mm).tif']));
    img_wht = im2double(img_wht);
    
    figure('units','normalized','outerposition',[0 0 1 1]);
    subplot(4,2,1)
    imshow(cb_prime)
    title('棋盘格原图');
    subplot(4,2,3)
    imshow(noise_prime)
    title('伯努利噪声原图');
    subplot(4,2,5)
    imshow(black_prime)
    title('黑板原图');
    subplot(4,2,7)
    imshow(white_prime)
    title('白板原图');
    subplot(4,2,2)
    imshow(img_cb)
    title('棋盘格拍摄图');
    subplot(4,2,4)
    imshow(img_noise)
    title('伯努利噪声拍摄图');
    subplot(4,2,6)
    imshow(img_blk)
    title('黑板拍摄图');
    subplot(4,2,8)
    imshow(img_wht)
    title('白板原图');
    pause(2);
    close;

    % 4.2 进行明暗校正
    u_img = img_blk + iprime_norm .* (img_wht - img_blk);
    figure()
    imshow(u_img)
    title('明暗校正后的潜在清晰图像')
    pause(1);
    close
    
    % 4.3 进行PSF估计
    this_shift = All_Position(i);
    R = magnify*abs(this_shift-focus_shift); % 模糊半径
    k_size = ceil(3*R);
    % 将k_size变成奇数
    if mod(k_size,2)==1
        k_size = k_size+2;
    else
        k_size = k_size+1;
    end
    k_size = max(min(psf_size,k_size),15);

    % Load the u_imgs and b_imgs for PSF estimation
    b_img = img_noise;
    

    u_img = u_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);
    b_img = b_img(rect(2):rect(2)+rect(4)-1,rect(1):rect(1)+rect(3)-1);

    k_init = zeros(k_size); k_init(ceil(k_size/2), ceil(k_size/2)) = 1;   

    tic
    psf = PSF_estimation(k_init,u_img,b_img,R);
    toc
    % 在填充前保存psf
    savearray{i} = psf;
    
    % 4.4 对PSF进行边缘零填充
    psf = padarray(psf,[(psf_size-k_size)/2,(psf_size-k_size)/2],'both');
    figure('units','normalized','outerposition',[0 0 1 1])
    imshow(psf,[]);
    filename = fullfile(imgpath,[num2str(wavelength),'nm.png']);
    exportgraphics(gcf,filename);
    pause(1);
    close;

    % 4.5 进行PSF去模糊测试结果
%     img = imread(fullfile('raw_data\20220117', '分辨率板' ,['分辨率板 ' num2str(wavelength) 'nm.tif']));
%     img = im2double(img);
%     img_deblur = deconvlucy(img,psf);
%     figure()
%     subplot(1,2,1)
%     imshow(img);
%     title('拍摄图');
%     subplot(1,2,2)
%     imshow(img_deblur);
%     title('去模糊图');
%     pause(1);
%     filename = fullfile(img1path,[num2str(wavelength),'nm.png']);
%     imwrite(img_deblur,filename);
%     close;
    
    % 4.6 保存数组(每隔10次)
    
    if mod(i,10)==0
        save(fullfile(savepath,'PSF_550nm.mat'),'savearray');
    end
end
